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, The estimator proposed recently by Delmas and Jourdain for waste-recycling Monte Carlo achieves 

variance reduction optimally with respect to a control variate that is evaluated directly using the 
simulation data. Here, the performance of this estimator is assessed numerically for free energy 
calculations in generic binary alloys and compared to those of other estimators taken from the 
literature. A systematic investigation with varying simulation parameters of a simplified system, 
the anti-ferromagnetic Ising model, is first carried out in the transmutation ensemble using path- 
sampling. We observe numerically that (i) the variance of the Delmas- Jourdain estimator is indeed 
i-S^ , reduced compared to that of other estimators; and that (ii) the resulting reduction is close to the 

^ ' maximal possible one, despite the inaccuracy in the estimated control variate. More extensive path- 

sampling simulations involving a FeCr alloy system described by a many-body potential additionally 
show that (iii) gradual transmutations accommodate the atomic frustrations, thus alleviating the 
numerical ergodicity issue present in numerous alloy systems and eventually enabling the determi- 
d ' nation of phase coexistence conditions. 
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In Monte Carlo simulations of multi-particle systems, thermodynamic quantities are traditionally estimated by an 
arithmetic mean taken over the Markov chain of states constructed by the simulation. As a result, the information 
pertaining to the trial states that have been generated and rejected by the Metropolis-Hasting algorithm is lost. 
^ ■ In contrast, waste-recycling Monte Carlo is a technique which includes this information in the estimate in order to 
1—1 1 improve its accuracy^ 

Waste-recycling estimators have been provided for quantum Monte Carlo schemes^ multi-proposal algorithms^ 
parallel tempering (with pair-replica exchanges^ or multiple-replica exchanges^), path-sampling^ and path-sampling 
, with multi-proposalj£~— In most implementations, waste-recycling was shown to be beneficial for estimating free 
f**- ' energies or potentials of mean force, in the sense that variance reduction was observed numerically. Cases of variance 
augmentation have however been reported ) 11 ' 12 showing that the additional recycled information is not always relevant. 
Nevertheless, Delmas and Jourdain demonstrated mathematically that variance reduction is indeed achieved in the 
\ asymptotic limit when the acceptance function is symmetric to the identity of the old and the proposed state, as is 
■ the case for the sampling algorithm named after Boltzmann, Glauber or Barker. The estimator including information 
from both the old and the proposed state is then the associated conditional expectation^ The authors also cast 
waste-recycling Monte Carlo into a general control variate problem structure and derived the optimal choice of the 
^ ' control variate in terms of asymptotic variance. 

The purpose of this article is to investigate the relevance of the optimal estimator for calculating free energy 
differences, a crucial task of molecular simulation^ The article is organized as follows. The multi-particle system is 
presented formally in Sec. [TT] along with a broader overview of the weighted path ensemble approach employed, while 
the more specific aspects concerning the simulated binary alloys are deferred until Sec. [V] Weighted path ensemble 
averages are derived in Sec. IIIII and the Monte Carlo algorithm, including the mono-proposal sampling scheme and 
various estimators, are described in Sec. IIVI In Sec. [V] the methodology is applied to the calculation of chemical 
potential differences in binary alloys. While the primary goal of the study is to compare the statistical variance of 
the Delmas- Jourdain estimator with that of other estimators, we also illustrate the possibilities of the methodology 
by determining coexistence conditions using the equal-area construction on the chemical potential surface or the 
common-tangent construction on the reconstructed Gibbs free energy surface. 
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II. EXTENDED SYSTEM AND DYNAMICS 



In this article we are interested in calculating the free energy difference for converting one particle of an alloy 
from one component type to the other, the systems and 1 discussed below. In order to achieve this, we will employ 
extended systems which allow us to convert the Hamiltonian of the system of one composition to the Hamiltonian of the 
system of the other composition via a switching parameter A, in a manner reminiscent of thermodynamic integration 
or nonequilibrium work calculations. In this section, we introduce much of the basic notation and terminology, and in 
the subsequent section, we discuss how the statistical physics of these paths generated while switching A may be used 
to construct an alternate work-weighted path ensemble which is particularly well-posed for transmutations between 
differing alloy compositions. 

The free energy of a multi-particle system of interest, labeled 1, can only be computed as a difference with respect 
to a reference thermodynamic state labeled 0. Let H\(x) and Hq(x) denote the Hamiltonians of the respective systems 
with x = (q,p) denoting the particle positions q and the particle momenta p composing phase space. We define an 
extended Hamiltonian H(x) = (1 — \)Ho(x) + XHi(x) for < A < 1 with respect to the extended state x — (\ x )- The 
switching parameter A is considered as an additional coordinate for the sake of generality, as the extended Hamiltonian 
can be given more general forms. The phase space associated with all extended states x sucn that A = a is denoted 
uj a , and the union u> = Uo<a<iWa defines the extended phase space. The probability densities in uj a and w are written 
respectively as 

Pa(x) = S x _ a ( X )exp\J3(G a -H a (x))] (1) 
p( X ) = exp[P(G-H( X ))} (2) 

where G a and G are the Gibbs free energy of the particle system and of the extended system, respectively, and the 
inverse temperature is j3. The pressure P and volume V are implicitly accounted for by the Hamiltonian. The delta 
measure 5\_ a (dx) = S\- a (x)dx is such that under Lebesgue integration J ip(X)5\_ a (dx) = <f( a ) f° r an Y test function 
ipJS. The Gibbs free energies G a and G are related to each other via the identity 

G a = G- /T 1 In / p( X )d X . (3) 



In the following, the switching parameter will continuously evolve between states and 1, but the two values of a 
of interest for calculating thermodynamic values will be and 1. The thermodynamical expectation of quantity <j) is 
expressed as 

(p a ,<f>) = J H x )Pa(x) d X- (4) 

Equation Q will be estimated using a waste-recycling Monte Carlo technique from Ref. [ill with appropriate 
modifications to reflect the recently developed Delmas-Jourdain optimal estimator. To improve ergodicity in the phase 
space as well as the accuracy of the estimates, the trial moves are trajectories generated within a path ensemble using 
Langevin dynamics at temperature and constant pressure P, yielded by coupling the particles to a thermostat 
and barostat. In addition, an external force /^ xt (x, A) depending on the extended state acts upon the additional 
variable A mechanically. In the most general set-up, A is equipped with a mass m\ and evolves according to 

m x X = f^ t -d x H(x,X), (5) 

assuming that A is not restricted to the [0,1] interval. A path z = {x*}o<t<r °f duration r consists of the sequence of 
the states Xo<*<r generated by the Langevin dynamics starting from xo- TTie conditional probability to generate path 
z knowing the initial state xo of the system is denoted by Pf w d(xo<t<r Ixo) and can be evaluated numerically from 

the normal random deviates used in the stochastic dynamics. The reverse of path z is denoted by z^ = I x\ \ 

I J 0<t<T 

and xl — {\ Qi ~P)r-t where (X,q,p) t = Xt- The conditional probability to generate path z backward knowing the 
final state Xr of the system is the probability to generate z^ forward from Xo 

Prov(XO<t<r|Xr) = Pfwd (xl< t <r IXo) • ( 6 ) 

and is also evaluable from the normal deviates required to generate z' . In the following, the labeling of the probabilities 
with fwd and rev is dropped and we rather write 



(Xo< t<ar \Xar )Pfwd(XaT<t<-r|Xa7-) 
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The two values of a of interest will be and 1. Hence, the nature of the probabilities is implied by the conditional 
presence of xo or Xt- 

The forward and reverse conditional probabilities above are two crucial quantities in waste-recycling Monte Carlo 
as their ratio enters the acceptance rule used both by the sampler fSec. HV A"| and the estimator fSec. UVBj) . The ratio 
is related to the heat Q(z) transferred from the thermostat and barostat into the system along z via the well-known 
expressio n 16 ' 17 

In presence of an external force acting upon A, the heat may be expressed as (see Eq. 26 and Eq. 27 in Ref. fToh : 

,A(t) 

Q(z)=H( Xt )-H(xo)- fr\Xt)X(t)dt. (8) 

J\(0) 

The integral of the product of the external force acting upon A by the displacement dX = Xdt corresponds to the 
mechanical work W{z) done on the extended system. As a result, the identity in Eq. (0 can be cast into the equivalent 
form 

P(z|Xt) exp [-PH(xt)] r aur , \i , n , 

—— — — - = exp [-pW(z)] . 9) 

P(z|xo)exp [-/3H(xo)\ 

that will prove convenient for waste-recycling Monte Carlo because it contains the Hamiltonians of the distribution(s) 
of interest. Note that the identities in Eq. ([7]) or Eq. © remain valid when the evolution equation Eq. (J5J) is coupled 
to a thermal bath via an Ornstein-Uhlenbeck process^ 

An autonomous scheduling of the additional coordinates based on the general evolution equation Eq. <j5j) enables 
the dynamics to explore all regions of interest while equilibrium information is retrieved using the ratio in Eq. Q 
by waste-recycling Monte Carlo. This approach is particularly relevant when free-energy is to be reconstructed in 
multiple dimensions, or when the additional coordinates are meta-variables ^ 18 ' 19 Meta-variables are associated with, 
for instance, a restraining harmonic spring acting upon a generalised or reaction coordinate and characterize the 
position of a pulling device added to the particle system. Although the additional device modifies the Hamiltonian of 
interest, its contribution to thermodynamic expectations can be canceled by taking the ensemble average directly in 
the extended ensemble^, defined as follows 

(p,tf>) = J cb(x) P (x)d X . (10) 

Meta-variables are most commonly employed when constructing the free energy along reaction coordinates via umbrella 
sampling, which in its usual implementation 14 entails correcting for the sampling bias rather than resorting to a 
marginal expectation like Eq. (fTU|) . 

Here, we are interested in thermodynamic states and 1 and need not consider Eq. (110[) . We further simplify 
Eq. (O by requiring that A varies at constant velocity from A(0) = to A(r) = 1, as advocated by Jarzynsk L 20 ' 21 The 
external force satisfying the resulting constraint A = in Eq. §5§ is such that /^ xt = d\H, hence the work done on 
the extended system is 



W{z) = / d x H{\ u x u )X t dt. (11) 



o 



III. WORK- WEIGHTED PATH ENSEMBLE 



In this section, we show how the extended systems of the previous section may be employed to construct the 
work-weighted path ensemble used in our simulations and sampled by waste-recycling Monte Carlo. 

In the following, the full trajectory space f2 encompasses all paths that connect the wo and u>\ phase subspace 
as the additional coordinate A varies monotonically from at t = to 1 at t = t. Besides, two parameters a and 
9 are employed which formally may adopt the full range of possible values from to 1. The intent of each of the 
parameters is distinct however. As previously mentioned, a is meant to indicate the thermodynamic states that we 
are interested in probing, and as such the meaningful values of a are strictly and 1. The other parameter 9 is a 
weighting factor supplied in constructing the path ensemble, in a sense indicating for each trajectory the contribution 
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of the two possible generating positions \o — (0, x$) and Xr = (T^r) to its associated probability. In the results, a 
full range of 8 values will be explored, but 9 = 0.5 is often advantageous, 5 allowing for a good overlap^ with the two 
work histograms that can be constructed when paths are generated either forward from the equilibrium distribution 
of system or backward from the equilibrium distribution of system 1. 

The generalized path ensemble has a weighted probability density Pg for the generating parameter 9 in O given 
asi^ 

Pg(z) = e« {P(z\xo) exp (xo)]} 1 ^ {P(z\Xr) cxp [-pH{ XT )]} (12) 

where the normalizing constant gg ensures that Pg is a probability distribution 

/ Pg{z)Vz = l. (13) 
Jn 

In the following discussion, the notation \ aT characterizes the state of path z at t = cut, allowing us to develop the 
very similar equations for the states of interest, xo an d Xn with more compact notation. For the two particular values 
a e {0, 1}, Eq. (jl3]l is equivalent to 



I e -P H (x) dx= i ( 14 ) 



which results from the simplification of Eq. (fl2|) into P a (z) = e 9a e ^ H ^ XaT ^P (z\xar) and from the normalization of 
the conditional probabilities 

/ P(z\ XaT )Vz = l (15) 

in each subspace Q(xar) consisting of paths going through Xar- As a result, the logarithm of Eq. ([T4|) with a E {0, 1} 
reads 

g a = - In f exp [-0H{ X )] d X = PG a , (16) 

indicating that these normalization constants go and g\ are indeed free energies for the two endpoint systems. 

The normalization of the conditional probabilities for a G {0, 1} also enables one to express the thermodynamic 
expectation of any quantity (f> with respect to the path density P Q rather than the state density p a (for a = or 1) 



( Pa ,4>) = I <l>(x)exp\P(G a -H(x))]dx 

S x (x a r)Hx)e^-^P a (z\x aT )Vz)d X 

! / 

= J 4>(x aT )P a (z)T)z 

= (Poi^a) 

where z — > 4> a (z) is the path functional such that 4> a {z) — 4>{x aT ). 

Since importance-sampling is achieved with respect to the distribution Pg with < 9 < 1, we should rather employ 
to the formal path-average 

(p a ,4>) = (P e ,0 Q P Q /P e ). (17) 

The probability density ratio in Eq. (|17p can be cast into the form 

P a /P e =e^ e -^ w+9 «- 9e 1 (18) 

obtained after substituting a for 9 in Eq. (fT2"|). dividing by Eq. (fTS"]) while leaving 9 intact, and further utilizing the 
connection between the conditional probabilities and W as defined in Eq. Furthermore, substituting e 9e ~ 9a for <p 
and 4>a m both sides of Eq. (fT7f yields 

P e ,e^ B -^ w ). (19) 
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The form of Eq. (|19[) enables one to express Eq. (|18[) without the unknown normalizing constants g a and gg, as follows 

P a /Pe - ef>V- a W w / (P e y (6 - a)w ) ■ (20) 

Inserting the probability ratio in Eq. ([20)1 into the path average of Eq. (fTT)) finally yields for a € {0, 1} 

<Pa,0 = (Pe^ a e^ w )/(Po,e^ w ). (21) 

The ratio involves two computationally tractable forms that can be estimated using Markov chain Monte Carlo 
methods. In particular, the free energy difference G\ — Gq = /3 (<?i — go) can be obtained from the relationship 

Gl - Go =-^^r xm LT ]) - (22) 

(P0,exp[/W]) V 1 

The effect of the bridging parameter 9 on the numerical performance is discussed in Refs. andl22l and further studied 
in Sec. [V] Note that other work-biased distributions than Pg and other ensemble averages have been proposed in 
the literatur o 23 ! 24 and are reviewed in Ref. l25l . Here, one recovers Jarzynski's identity by choosing the importance 
function Pg to be Po 

Gi - G = -/T 1 In (P , exp \-pW\) . (23) 

Similarly, Crooks's probability ratio^l 

P(z\Xt)Pi(Xt) 



P(z\xo)Po(Xo) 



exp\fi(G 1 -G -W{z))]. (24) 



is recovered from Eq. (|20|) after setting 9 = and a = 1 and then substituting Ppo and Ppi for Po and Pi respec- 
tively Equation (1241) and Equation ([§]) are two mathematically equivalent identities referring to two distinct physical 
frameworks. In Eq. (|24[) . W(z) is to be interpreted as a nonequilibrium thermodynamic work, defined as the integral 
of the energies 8E(X; x) = H(X + SX, x) — H(X, x) transferred to the system as accounted for by small changes SX in 
the external generalized mechanical constraint A i 26 > 27 Constraining the system via the control parameter A amounts 
to specifying the equilibrium distribution p\. 

Within this framework, the identity in Eq. (|24l) entails a free energy difference. In contrast, with the physical 
prescription of Sec. HU the external mechanical force /^ xt acts transiently on A for a duration r and is then switched 
off. As a result, the system returns to its equilibrium distribution p(x) = ex P [G — /3H (x)]- In this alternative 
framework) 2 ^ the reverse and forward occurrence probabilities of path z in the ratio in Eq. © rather corresponds to 
a Bochkov-Kuzovlev identity . 28 ' 29 

Neither of the two above frameworks holds with respect to the weighted path ensembles with distributions P# when 
< 9 < 1 because, in this situation, no Hamiltonian-based equilibrium distribution can be defined for the state 
ensemble. The path formalism is merely introduced to facilitate the exploration of the phase space as will be shown 
in Sec. El We now discuss how to sample Pg and to estimate (Pg, ■) via a Markov chain Monte Carlo algorithm. 



IV. WASTE-RECYCLING MONTE CARLO WITH MARKOV WEBS 



A. Samplers 

In practice, the distribution Pg is approximated by a Markov chain constructed by importance sampling. Any 
sampler consists of two steps: (i) starting from a given path Zk, a trial path Zk is generated from a probability 
distribution q; (ii) the trial path is accepted and added to the Markov chain Zk+\ = Zk with an adequate probability 
p, otherwise Zfe+i = Zk- So as to ensure the convergence of the Markov chain towards the correct distribution, the 
transition probability matrix P must satisfy the following detailed balance equation 

P{z s ,z)P e (z) = P(z,z s )Pg{z s ), (25) 

where P(z s , z) is the probability to transit from path z to path z s , and vice versa for P(z, z s ). 

To formalize the sampler, we write the probability to construct the proposal z from z as g({z,z}|z) and the 
acceptance probability of the proposal as p(z,{z,z}), making the set of the proposed path and the original path 
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explicit in each case. The rejection probability is 1 — p(z, {z, z}) and amounts to transiting to z. The transition 
probability from z to z s is therefore^ 3 - 

P(z s , z) = J2 [6z(z s )p& S}) + S z (z s ) (1 - p(z, {z, z}))] q({z, z}\z), (26) 

z 

where 5 is the delta distribution and we allow for the possibility that z s is either the old path z or the proposed path 
z. We write g({z,z}|z) to express the backward probability enabling one to construct the current path z from the 
proposal z. Then, an acceptance probability that can be constructed from Eq. (|25l) and Eq. (|2"6"|) is 



p M (z, {z, z}) = min 1, ' ' P . (27) 
V Q({z,z}\z)P e (z) J 

This rule, proposed by Metropolis, is traditionally used because it maximizes the probability of acceptance. In 
particular, it is always larger than the acceptance probability of the Barker algorithm given by 

q{{z, z}\z)P e {z) + q({z, z}\z)Pg{z) 

This acceptance probability is symmetric so that p B {z\{z 1 z]) — p B (z\{z, z}), in contrast to the Metropolis rule p M . 
Moreover, p B appears to be a posterior likelihood associated with a marginal probability that is actually sampled by 
the algorithm (see Ref. [ll| and Appendix [5} ■ 

We consider that the sampler is such that the generated trial paths z share a common state Xar with the current 
path z where a £ {0,1}. The constructed pairs of such paths are denoted by {z, z} a and are called webs. The 
proposal functions are given by the conditional probabilities 

q({z,z} a \z) = r a P(z\ XaT ). (29) 

The factor r a relates to the probability to choose a = or a = 1. We similarly define q({z, z} a \z) = r a P(z\x aT ) 
for the reverse move, with the couple (ro,ri) in practice chosen to be alternately (0, 1) and (1,0). Note that another 
possible choice consisting of setting r a = \ both for a = and a = 1 appears less efficient numerically.— 
The Barker acceptance rule is therefore (a G {0, 1}) 

v B (~z\\z ~z\ a ) = exp[/3(q-g)W^)] 

P K [X , i 1 exp[l3(a-9)W(z)}+exp[p(a-e)W(z)Y [ ' 

Acceptance ratios for < a < 1 are derived in the literatur o 22 ! 30 but are not implemented here. 

Once the sampling process has been completed, the quantities of interest can be estimated from the constructed 
chain of webs {z k , h} 1<k<n - 



B. Estimators 



The standard estimator for the expectation (Pg, /) of a path functional / from a Markov chain of paths {-Zfe}i</ C <„ 
distributed under Pg is defined as 

1 - 

Uf) = (si) 

fc=i 

which is, in essence, the straightforward summation of states in the chain. In applications of Sec.|Vl the path functional 
/ will be chosen to be / = e^ e ~ a ) w . 

Waste-recycling estimators are based on conditional expectations 13 and include information from trial moves using 
the acceptance probability.— £ In the path-ensemble, the waste-recycling estimator reads 

n 

where the nature of the webs indicating a shared state at either a = or a = 1 is omitted for brevity and {z kl z k }i< k < n 
is the web chain. 
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Following Delmas and Jourdain^ waste-recycling is reformulated as a control variate problem consisting of search- 
ing the value & gl minimizing the asymptotic variances <Jg(f, bf) 2 of the estimators J^if) defined by 

J b n (f) = (l-b)I n (f) + bI™(f). (33) 

Note that with the present notation, we have J° = /„ and J\ = for the standard and waste-recycling estimators 
with respective asymptotic variances crg(f, 0) 2 and <Jg(f,f) 2 . Delmas and Jourdain proved that the function b — > 
a e{f,bf) 2 — <rg(f,0) 2 is a convex quadratic form whose coefficients are easily evaluable when sampling is performed 
with the Barker algorithm. To make this remarkable feature of ag(f,bf) 2 clear, we shall write Kg [1/3(^1,^0)] for the 
conditional expectation of function ip when zq is distributed under the invariant measure Pg and Z\ is a random 
deviate constructed from zq using the sampler whose transition probability is P(-, zq) given in Eq. (|26p . 
The obtained quadratic formic 

Mf, bf) 2 - Mf, 0) 2 = \k 6 [(/ (zi) - f (z )f] xb 2 - 2{(P g , f 2 ) - (PgJ) 2 } x b (34) 

is valid for Barker sampling only. The coefficient of the first order term in b is the variance of the observable. 
The coefficient of the second order term is a conditional variance since Kg[f(zi) — /(zq)} = (we have (Pg,f) = 
Kg [f{zo)] = Kg [/(zi)]). Since both variances are strictly positive for non-constant / on {z/ 'Pg(z) > 0}, the quadratic 
form Eq. (|34p has a minimum occurring at the strictly positive value 

b* = < p ;./ a >-< p */>V (35 ) 



\K e 



(/(^i)-/M) 2 



It is moreover proven^ 3 - that b* > 1 when <r(f, J) 2 > 0. To simplify the notation, the dependence of b* on 6 and / is 
omitted. 

This implies that the value b* that minimizes the variance can be estimated from a single run in spite of the fact 
the variance <7g(f, 0) 2 of the standard estimator can only be evaluated from a sample of estimates I n (f) taken in the 
limit of large n owing to the central limit theorem^ The optimal parameter b* defined by Eq. (|3"5)) can be estimated 
by 

X n{b ) = I n r ~ ■ 36) 

s;E*=i(/(**+i) ~f( z k)) 2 



from the Markov chain. An alternative procedure to estimate the optimal parameter Eq. (|35|) consists of in- 
cluding information about trial paths Zk- In particular, Kg[<j)(z-L, Zq)] can be estimated via waste recycling by 

~12k=iP B (zk\{zk, Zk})4>(zk, z k ) + p(zk\{zk,zk})(f>{zk,zk). With <j>(zi,Zo) = (f(zi) - f(z )) 2 , the optimal parame- 
ter b* defined by Eq. ([33)) can be estimated using 

I WR (b*) = (37) 

SI Hk=iP( S k\{zk, Zk}){f{z k ) ~ f{z k )) 2 ' 

Note that the waste-recycling estimator is not optimal at b* in terms of asymptotic variance when sampling is 
achieved with the Metropolis algorithm. The question whether one should resort to the optimal estimator and the 
suboptimal Barker sampler or to the optimal Metropolis sampler and a suboptimal estimator has not been answered 
theoretically and will be addressed numerically here. 



V. APPLICATIONS 



We now apply this waste recycling Monte Carlo on Markov webs to our system of interest, binary alloys of varying 
composition. Initially, we shall make clear how the general notation of the preceding sections is connected to this 
problem of interest. Subsequently, we delve in to the numerical applications. The performance of the optimal estimator 
is first assessed by comparing with that of traditional estimators in a generic but realistic binary system with A and B 
atoms interacting on a rigid lattice in Sec. IV El The possibilities of the methodology are then illustrated by performing 
off-lattice simulations of a model FeCr alloy in Sec. IV CI 
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A. Binary alloys 



The computational cell contains N atoms. In the reference state at t — 0, we have N = Na + Nb or N = Np c + Nci 
atoms where Nx refers to the atom number of type X. In the target state at t = r, the cells still contains N atoms, but 
an A or Fe atom has been transmuted into a B or Cr atom, respectively. Paths thus consist of artificially switching 
the potential energy of a selected atom. The switching will be performed instantaneously in the AB alloy fSec. IV Bl . 
or gradually using constant-pressure Langevin dynamics^ in the FeCr alloy fSec. IV Cj) . 

The infinitely- fast transmutations of Sec. IV Bl are carried out without changing the atomic masses. Hence, the ratio 
of the reverse conditional probability, P(z|x r ) = (Nb + to the forward one, P(z|xo) = (N — Ab)^ 1 , relates to 

the exponential of the ideal chemical potential difference 

Plugging Eq. (|38|) into Eq. ([9]) leads to an identity relating the work W carried out on the system for transmuting an 
A atom into a B atom to the Hamiltonian variation 

W(z) = H( X r) - H(xo) - A M id (39) 

where Xt only differs from xo by the transmuted atom. When transmuting a B into an A atom, the quantity —W[z) 
must be considered instead. Besides, A^ ld acts as a heat transfered from a reservoir of A and B atoms into the 
system. 

In contrast, the transmutations of Sec. IV CI are performed gradually using constant-pressure Langevin dynamics^ 
and linear Hamiltonian switching. As a result, the probability ratio Eq. (|38|) becomes 

^4=exp[/3A M id + /3Q(^)] (40) 

where Q(z) is the heat transferred from the thermostat and barostat to the particle system^ in addition to the heat 
A/x ld transfered from the atomic reservoir. We thus have 

W(z) = H( Xt ) - H( Xo ) - A^ id ~ Q(z) (41) 

In both set-ups, sampling forward and backward transmutations is sufficient to explore the phase spaces of the alloy. 
As illustrated in Fig. [TJ accepting several paths amounts to exchanging the allocation of atoms on the underlying 
lattice of the reference and target systems. 

Simulations aim at extracting differences of chemical potentials A^i between both species as a function of the alloy 
composition. We have either A/i = /j,b — MA or A/i = /io — //Fe depending on the involved alloy system. Both 
quantities indeed correspond to the difference of Gibbs free energy between the target and reference states 

A^ = /?- 1 ( 3l - 5o ). (42) 

The generic estimator will be 

jb ( j3(e-i)w\ 

J n 6 (A/x) = -/3- 1 ln j bn{efieW) } ( 43 ) 

The numerator and denominator in the logarithm are estimates of exp [gg — g{\ and exp [gg — go] obtained simultane- 
ously by sampling the probability measure Pg and appling the generalized estimator in Eq. (I33|) . 

Let c = Nb/N or c = Ncr/N denote the alloy concentration in B or Cr. A change in the monotonicity of the 
function c — > A/i(c) is a signature of phase coexistence. The conditions of phase equilibria can then be determined via 
the equal-area construction with respect to A/i(c), or equivalently, via the common-tangent construction with respect 
to the Gibbs free energy 

G(c) = f An{c')dc'. (44) 
Jo 



measured per atom. 
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FIG. 1: Schematic free energy landscape in the extended space. The alloy undergoes forward and reverse atomic transmutations: 
white curves linking free energy minima symbolize the sampled transmutations. Path-sampling enables the simulation to explore 
the phase space of both the reference and target systems. 

B. AB system with constant pair interactions 

In our simplified binary alloy lattice model, the rigid lattice is body centered cubic. Interaction energies are taken 
as pair interactions exy between nearest-neighbor sites, where X and Y may equal A or B. The ordering enthalpy 
e = f-AA + £bb — 2eab plays a key role as it entirely determines the thermodynamics of the system^ Without loss of 
generality, we can choose €aa — £ ab — which leads to an Ising-type internal energy 

H( X ) = e J2 VW B (45) 
(i,j)eE 

where the summation runs over the ensemble E of nearest neighbor pairs only and vj^ is 1 or depending on whether 
site i is occupied by a B atom or not. The cell contains N — 2 12 sites with JVb atoms of type B at t = 0. 

We set e = —30 meV. With a negative ordering energy, the system exhibits a miscibility gap below which the solid 
solution decomposes into A-rich and B-rich phases. The unmixing transition is of first-order except for the A0.5B0.5 
composition where it is second-order— and where the critical temperature is T c m —e/(k x 0.62)M- 

We first check the Delmas-Jourdain prediction that the asymptotic variance a(f, bf) 2 of J^if) is minimal at b = b* 
when the Barker sampler is used. Since lim nvar„(/, bf) — <j(/, bf) 2 where var„(/, bf) is the statistical variance of 

71— > + OG 

J„(f), we evaluate var n (/, bf) as a function of b for large enough n and confirm whether its minimum occurs close to 
estimated values of b* . Note that b — > var„(/, bf) is a quadratic form regardless of the value of n. 

We set temperature to 348 K, the simulation parameter 9 to 1/2, and B concentration to 50at.%. The two quantities 
exp [<7i/2 — go] an d exp [gi/2 — <7i] are equal to each other for this particular set-up, and have been evaluated via 
ensemble average Eq. ([TTJ)) and the estimator J*(f) given in Eq. K33J) with / = e ± ^ w and < b < 20. The statistical 
variances var n (/, bf) have been computed from m = 10 7 estimates generated using distinct random seeds. Each 
estimate consists of n = 2 ■ 10 3 transmutations, from A into B and from B into A alternatively. For each estimation, 
we also record the estimated value of b* using estimator X n (b*) given in Eq. (HU and estimator I™ R (b*) given m 
Eq. (|3"T|) . and additionally construct their histograms: (h°(b — &*)} and (h 1 (fe — b*)) respectively denote the probabilities 
that estimates l n (b*) or l^ R (b*) are equal to b. 

As shown in Fig. f2J the quadratic form var„(/, bf) is indeed minimum at the value corresponding to the best b*- 
estimate indicated by the vertical double-dotted segment, obtained by combining the mn available data from either 
3~nm(b*) or T™^(b*). The minimum is also very close to the horizontal line labeled var n (/, /&*) corresponding to a 
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FIG. 2: (left ordinate) var„(/, 6/) is the statistical variance of J^(e l3w ^ 2 ) given in Eq. (|33[) and is evaluated from 10 7 estimates. 
Note that var n (/, 0) and var„(/, /) corresponds to the variances ofi„(e-W 2 ) and/r(e- WJ ). 

var„(/, 6*) is defined in the 

text; (right ordinate) (h°(6* — b)) and (h 1 (6* — b)) are the histograms of optimal paramameter b* estimated without or with 
waste-recycling. 



variance obtained from the m estimates </«"(/) after substituting the corresponding I n {b*) for &* in each estimate. 
Substituting 2^ R (b*) for 6* further decreases the variance by 0.24%, an amount not visible on the graph. A more 
noticeable benefit to employing waste-recycling for estimating b* can be seen from the two histograms of b* displayed 
in Fig. [2j Histogram (h 1 (6 — b*)} is slightly narrower than histogram (h°(b — b*)) obtained without waste-recycling. 

We also observe in Fig. |2] that the optimal estimator exhibits a smaller statistical variance than that of the residence- 
weight estimator K n (derived in Appendix IA")). although the available theory docs not predict the extent to which 
the observed tendency is general. We also don't know the value 6* of the control parameter that minimizes the total 
variance in the estimates of A/x from Eq. (|4"3"1) . Qualitatively, the symmetric setting 8 = 1/2 is often used&£ ' 22 i 30 
because the work distribution associated with P1/2 exhibits sufficient overlaps with those associated with both Po 
and Pi, ensuring a fast convergence of the involved exponential averages^ 

We carried out a series of simulations with varying the value of 8 in the range [0,1] so as to locate the optimal 
value 8* with respect to the evaluation of A/x using estimator Eq. (|43p . at the asymmetric composition c = 10at.% 
B. The statistical variance is obtained using m = 10 7 independent estimates but with n = 2 • 10 4 transmutations 
in each estimate (10 times more than previously). The results in Fig. [3J demonstrate that the observed hierarchy in 
the estimator performance is preserved for any 8 6 [0,1]. The optimal estimators achieve the best performance 
whether waste- recycling is used (h 1 ) or not (h°) for estimating b* . Interestingly, using the optimal estimator shifts 
the minimum of the variance to a value of 8 close to 0.5 from a value about 0.7 with the standard estimator. Given 
that b* is about 3 when 8 = 1/2 close to the minimum of the variance, the "standard" contribution involving the 
Markov chain of paths accepted in the sampled distribution Pg is assigned with a weighting factor 1 — 6* close to the 
negative value of —2. The "waste-recycling" contribution involves both the accepted and trial paths and is assigned 
with a weighting factor b* that is positive and larger than 1 — b* in absolute value. Hence, compared to the overall 
contribution of the accepted paths, trial transmutions generated from the forward and backward path distributions 
contribute to the free-energy estimates quite substantially. 

We have then post-processed the collected data (the 10 7 estimates) by evaluating the statistical variances var ra (/, bf) 
of estimator J*(/) (given in Eq. ([33]) with / = ^l^-^W ) a e {q, 1} and n = 2 • 10 4 ) as a function of b 6 [0, 20]. 
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FIG. 3: (a) estimator variances; (b) acceptance rates; (c) optimal paramater; (d) reduced deviation of 6* as a function 9 
[var n (b*) is the variance of the estimates of &*]. Simulations are performed at T — 348 K. 



The mimima of the quadratic forms in b are plotted in panel (c) of Fig. [3] as a function of 9 (curve labeled 'post- 
processing'). As expected theoretically, they perfectly coincide with the online estimations of b* without waste 
recycling as in Eq. (|36|) or with waste recycling as in Eq. (|37|) . which are also plotted for comparison. Panel (d) of 
Fig. [3J shows the higher accuracy obtained for the waste-recycling estimates via Eq. ([57|) . 

Additional simulations are carried out at temperature /3 _1 = 174 K with varying 9 or at f3~ 1 — 348 K with varying 
concentration (and imposing 6 = 1/2). Results are displayed in Fig. 2] and Fig. O respectively. The same qualitative 
trends are observed. Nonetheless, for the low temperature simulation, much smaller acceptance rates are measured. 
A small acceptance rate decreases the denominator in both Eq. (j3l)|) and Eq. (|37|) explaining the higher values of 
b* and the difficulties in estimating the variance when 9 < 0.25 where a few negative estimates of e 9e ~ 9c " prevented 
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FIG. 4: Same as in Fig. E] but at T = 17 A K. 

evaluation of the logarithm in Eq. (j43|) . 

We have carried out additional Monte Carlo simulations using the Metropolis sampler with 9 = 0.5 and /3" 1 = 
348 K and with varying B concentration in order to compute the minima of the statistical variances of estimator 
J^(e /3( - 1 / 2_Q! ' IM/ ') as a function of concentration and for a € {0,1}. Similarly, 10 7 estimates have been used. The 
minima 6 mm are plotted as a function of concentration in Fig. [5] (curve labeled 'post-processing') and do not coincide 
with the online estimation in Eq. (|3"o| or Eq. (|3"T)) . as expected for the Metropolis sampler. 

Because evaluating statistical variances requires a considerable computational cost, accurately determining b mm by 
post-processing from the minimum of the function b — > var„ (/,&/) is more difficult than accurately estimating /, 
the quantity of interest. Post-processing the data obtained using the Metropolis sampler is thus not suitable to the 
simulation of realistic systems for which accuracy is an issue. 
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FIG. 5: (a) disagreement between b* and b mm obtained from the true minimum after post-processing with the Metropolis 
sampler; (b) agreement between the same quantities with the Barker sampler; (c) statistical variances of the 6* estimated with 
the Barker sampler. 

The question of whether the combination of the optimal Metropolis sampler and the suboptimal estimators /„ or 
R may achieve better efficiency than the combination of the suboptimal Barker sampler and the optimal estimator 
is not answered by theoryA2 Therefore, we now address this numerically. The statistical variances of the A/i- 
estimates obtained using the biased estimator of Eq. (H31 have been calculated and plotted as a function of B 
concentration in Fig. [5] We observe that it is always more efficient to implement the optimized estimator with 
the suboptimal sampler than estimator J\ with Metropolis sampling, or even the residence-weight estimator K n 
with Metropolis sampling. Note that is still a valid and unbiased estimator when combined with the Metropolis 
sampler (but is not optimal anymore). Here, the combination was found slightly more efficient, with variance further 
decreased by 4.7% for Aq ^Bq ^ compared to Barker sampling and i b n . However, variance reduction is not guaranteed 
mathematically and cases of variance augmentations have been reported when implementing waste-recycling estimators 
with a Metropolis sampleri 11 ^ 3 

Hence, the combination of Metropolis samplers with estimation of b is not recommended, is not presented in the 
plots, and is not further considered. 
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FIG. 6: (a,b) statistical variances of estimated differences of excess chemical potentials [displayed in (c,d)] as a function of 
concentration. Used estimators are indicated in the legends. The symmetry of this Ising system with respect to Cb = 50% 
explains the equality between chemical potentials (A/i = 0) at the solubility limit C s . 



C. Fe-Cr system with EAM interactions 



We now test our estimators on a more difficult model whose inter-atomic potentials are based on the embedded atom 
method (EAM) . We have implemented the EAM potentials of Olsson et aL— developed to model the a and a' phases 
of the FeCr binary system. While this EAM potential correctly reproduces the BCC structure of Iron and Chromium, 
its ground state (found among 99 candidate BCC structures based on the theory of the convex polyhedron in the 
correlation functions spaced) is an intermetallic structure with positive ordering enthalpy at c—50 at.% Cr^ So as to 
avoid the formation of additional phase stability fields, the intermetallic structure was excluded from the construction 
of the complete phase diagram using the Alloy-Theoretic Automated Toolkit (ATAT) package in Ref. |35|. The ATAT 
approach involves, firstly, constructing effective cluster interaction energies on a rigid lattice based on a few relevant 
structures previously relaxed using the EAM potential and, secondly, performing Monte Carlo simulations on the rigid 
lattice. Note however that the expected miscibility gap of FeCr alloy was observed in Monte Carlo simulations carried 
out on a rigid lattice with the interaction energies directly deduced from the EAM potential^ Here, the rigid-lattice 
assumption is entirely released in the construction of the phase diagram. 

We perform path-sampling simulations to compute the chemical potential difference with varying concentration 
and temperature. Transmutations are now performed gradually in 10 2 steps with linear Hamiltonian switching and 
constant-pressure Langevin dynamics 2 - in a computational cell containing 432 atoms. Equilibration proceeds in 20 
transmutations per atom, starting from a random distribution of the atoms on 6 3 unit cells of the BCC structure. 
Simulations have been carried out with temperature ranging from 300 to 1700 K in step of 25 K. Examination of 
the simulated microstructures show evidende of phase separation at low enough temperatures. Figure [7] displays the 
snapshot of a typical phase separated microstructure for alloy Feo.sCro.2 at 300 K. 

Figure [3] displays the statistical variances obtained for the various estimators obtained with n = 400 transmutations 
(per estimate) and m = 200 estimates at T=500 K (T — (fc/3) -1 ). The aforementioned hierarchy still holds in the 
present case for all concentrations. The Delmas-Jourdain estimator still achieves better efficiency than the residence- 
weight estimator. Because the amount of transmutations per estimate is much smaller than previously, it is more 
relevant to estimate b* using waste- recycling (h 1 ) than without (h°). The former variant further decreases the variance 
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FIG. 7: Snapshot of a simulated microstructure containing 20% at. Cr obtained at T — 300K. Fe and Cr atoms are displayed 
in blue and gray, respectively. 

associated with the estimation of b* by 20 % on average over the latter one. Averaging over all temperatures, estimator 
-h 1 decreases the variance by a factor of 2 to 5/2 compared to J® and by 20 to 30 % compared to with h°. 
The A/x-values displayed in Fig. correspond to a single estimate obtained from the 8 • 10 4 = m x n transmutations 
used to evaluate the statistical variance. Figures [8}d and[8j; illustrate the equal- area or common-tangent constructions 
of Maxwell and corroborate the occurrence of phase coexistence. 

The reconstructed Gibbs free energy surface G is less fluctuating because it is based on a self-averaging integral 
Eq. (EH)) . The potential quantity A/u* (T) in Fig. [5J: is such that the occurrence probabilities of the Fe-rich and Cr-rich 
phases are equal. The occurrence probability is P*(c,T) = Z^e^^*" ^ where Z = J* e N ^ c ' A ^- G( ~ c '^dc' is the 
semigrand canonical partition function* 4 associated with the finite computational cell. The potential function A/z*(T) 
exactly coincides with the slope of the common-tangent at coexistence and in the thermodynamic limit only. Within 
finite computational cells, Afi*(T) is defined even when there is no common tangent, is always easy to determine 
and fluctuates less than the slope of the common-tangent when there exists one. The transformed Gibbs free energy 
surfaces G(c) — cA/i* are displayed in Fig. |H1 The contour plot in the temperature-concentration plane of the bottom 
panel clearly visualizes the miscibility gap in the composition range going approximately from Feo.iCro.g to Feo.gCro.i. 

The coexistence lines (solubility limits) associated with the gap end around 500 K both in the iron-rich or chromium- 
rich sides and are indicated by the solid white curve in Fig. Above 525 K, we observe that free energy profiles 
becomes lower for intermediate composition which indicates the presence of a stability field for these intermediate 
compositions and of two immiscibility fields for the Fe-rich and Cr-rich compositions. Solubility limits on the Fe-rich 
and Cr-rich sides associated with the two immiscibility fields are indicated by the dashed white curves in Fig. [9Jd. 
However, examination of snapshots of the simulated systems within the expected stability field still shows unmixed mi- 
crostructures with a tendency to phase separation decreasing with increasing temperature from 525 K to 1700 K. This 
is attributed to the fact that the correlation length diverges at the critical temperature T c where the thermodynamic 
transition is second order. Given that the computational cell contains only 2 x 6 3 atoms and that atomic interactions 
up to 7th nearest neighbors were shown to play an important role^ strong finite-size effects are expected. The low 
temperature of 500 K measured for the present closure of the miscibility gap should therefore not be interpreted as an 
estimate of the critical temperature T c whose experimental value is expected to lie around 900 A finite-size scal- 
ing analysis 3 should therefore be carried out, which entails performing simulations with much larger computational 
cells. To achieve this task, the overall computation time has to be reduced considerably, possibly by evaluating the 
interatomic potential and forces on parallel computer architectures/ 3 ^ While the preliminary results presented here 
show that a direct and accurate construction of the equilibrium phase diagram of FeCr alloy is achievable in principle, 
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FIG. 8: (a) statistical variances of estimated differences of excess chemical potentials [displayed in (b)];(c) solubility limits 
displayed in the concentration-temperature plane as obtained from Maxwell equal-area construction at T = 500 K. Used 
estimators are indicated in the legend. 

they also emphasize the need for more extensive free energy calculations in this system. 

VI. CONCLUDING REMARKS 

In this study, the optimal estimator proposed by Delmas and Jourdain for waste-recycling Monte Carlo has been 
assessed numerically. As a testbed, we implemented the mono-proposal Barker and Metropolis-Hasting algorithms 
in transmutation ensembles of two alloy systems simulated via the weighted work path ensemble. We simultaneously 
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FIG. 9: Top are the Gibbs free energies transformed for better visualization (a = 10 -3 and b = hxZ/N) as a function of 
concentration for various temperatures; bright regions on the bottom contour plot correspond to stability fields, the Fe and Cr 
solubility limits of the Fe-rich and Cr-rich phases are schematized by the two white lines. 



estimated the free energy differences from the nonequilibrium works measured along the transmutations. We find that 
the estimator indeed achieves variance reduction compared to the other Monte Carlo estimators that are compatible 
with the present sampling approach. Furthermore, the maximal reduction of the statistical variance that is predicted 
by the theory is attained for relatively short simulations (2 • 10 3 sampled paths). 

Two other important methodologies have been proposed in the literature for optimally extracting free energy 
differences from nonequilibrium work measurements. Because they both require implementing specific sampling 
schemes not suited to the present alloy systems, we did not make any numerical comparisons. The differences in the 
implementation can nevertheless be emphasized. 
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The first methodology, proposed by Crooks^ and extended in Refs. [3^ and 0, entails sampling the reference 
and target thermodynamic states independently while performing bidirectional nonequilibrium work measurements 
starting from the sampled states. One then resorts to the optimal estimator of Bennett's acceptance ratio method.— A 
post-processing procedure constructs the optimal estimator from Bayesian inference, neglecting correlations between 
data4^ As the procedure minimizes the variances obtained from a finite sample with respect to some control parameter, 
the resulting estimator exhibits a statistical bias which becomes negligible only in the limit of large samples, whereas 
the optimal estimator of Delmas and Jourdain may be stably constructed when not in this limit and is always unbiased. 

One limitation of this alternative method relates to the lack of diagnostic tool to check the overlap between the 
forward and reverse switching processes. Insufficient overlap may result in poorly converged estimates. Another 
limitation lies in excluding the possibility of swapping between the reference and target states. For instance, when 
the target state corresponds to a low temperature structure that traps the system in a few basins of attraction, a 
single simulation of the target system may not fully sample all of these basins42 Swapping between the reference and 
target state as employed in this paper can avoid such non-ergodic sampling. 

The second methodolog y 15 ' 44 entails constructing an optimal biasing distribution rather than optimally post- 
processing the data obtained in various thermodynamic states. The optimal biasing distribution is derived analytically 
by minimizing an asymptotic statistical variance given by the central limit theorem^ 5 - This approach possesses a num- 
ber of limitations which we argue that waste recycling on work-weighted transmutation ensembles does not. First of 
all, the central limit theorem assumes independent and identically distributed random variables while sampled data 
are inevitably correlated; the denominator of the Delmas- Jourdain control parameter b accounts for these correla- 
tions. Secondly, the rejected trial paths are discarded from the optimization procedure, potentially losing valuable 
information. Thirdly, the shape of the optimal distribution crucially depends on the free energy difference that one 
wishes to compute and cannot be analytically estimated a priori. And finally, the optimal bias distribution con- 
structed in this methodology is bimodal with respect to the work with its two peaks corresponding to the forward 
and reverse work distributions and with a density depletion in the overlapping region of the two latter distributions. 
This may result in highly correlated data and poor statistics compared to our work-weighted transmutation ensemble 
simulations that sample a unimodal distribution exhibiting good overlapping properties in term of work with respect 
to the distributions obtained with both forward and reverse switching. Overall, we find that the Delmas- Jourdain 
estimator is unbiased and optimal in terms of asymptotic variance with respect to a simple control variate b* that 
can be accurately estimated from the correlations in the collected data. One substantial benefit to employing this 
optimal estimator is the flexibility in the choice of the sampler, allowing us to couple waste recycling Monte Carlo 
with a work-weighted transmutation path ensemble ideally formulated to study the phase coexistence of binary alloys. 
However, in other scenarios, another Monte Carlo sampling scheme might be better posed, yet the optimal estimator 
of Delmas and Jourdain would still be applied quite similarly. Since information associated with trial switching pro- 
cesses is optimally included in the free-energy estimates for any sampler, this one is to be chosen so as to facilitate 
the phase space exploration. 

Nonetheless, our approach employing the optimal estimator of Delmas and Jourdain shares a qualitative feature 
with the methodology based on the optimal biasing distribution. 44 Indeed, forward and reverse trial paths, which 
contribute substantially to the free energy estimate given by the Delmas- Jourdain estimator for the measured optimal 
values of the control variate b*, concommitantly cover the two peaks of the optimal bias distribution. 

Concerning the investigated examples of binary alloys, we point out that achieving numerical ergodicity in the 
reference and target thermodynamic states entails exchanging the allocation of atoms of distinct types on the under- 
lying lattice (which amounts to performing pairs of simultaneous transmutations with opposite direction) with a high 
enough frequency. With path-sampling, exchanges of atom allocation are automatically achieved when trial paths are 
successively accepted within the transmutation ensemble (as illustrated in Fig. [1]). It turns out that the phase space 
exploration is considerably facilitated. Resorting to such a path-sampling scheme was found particularly advantageous 
in the FeCr alloy system which presents a large magnetic misfit. Direct exchange moves sampled using a standard 
scheme in the reference and target system would have been extremely infrequent and thus computationally expensive 
compared to the gradual transmutations considered here. The approach might work as well in the numerous alloy 
systems exhibiting large atomic misfits. 

Further assessments of the estimator should consider more general scheduling for the nonequilibrium dynamics 
rather than simply Eq. ([5]) and a multi-proposal sampling scheme. Optimal waste-recycling estimators can possibly 
be derived with multiple proposals based on the same control variate problem structured Work in this direction is 
in progress ^ 



19 



Acknowledgments 

We warmly thank Benjamin Jourdain for suggesting that we test the optimal estimator and for advising us. Stim- 
ulating discussions with Berend Smit, Christoph Dellago and Peter Bolhuis are also gratefully acknowledged. 



Appendix A: Residence-weight estimator 



In the context of path-sampling, waste-recycling estimators may also be derived from statistical inference. The 
algorithm indeed samples a marginal probability distribution associated with the Bayesian prior Pg. As a result, 
information relative to unselected trajectories or to states contained in the trajectories themselves may be inferred 
online from a Bayesian posterior likelihood without post-processing. The marginal probability density of the sampled 
web {z, z} is 



®e ({z, z}) = I [q{{z, z} \z)P e (z) + q({z, z} \z)V {z)\ 



(Al) 



recalling that q is the conditional probability to construct {z, z} from either z or z. The factor 1/2 acts as a normalizing 
constant and relates to the position of the prior z in {z, z}. P$(z) an d P#(<?) are the prior probabilities of {z, z}. 

With the Boltzmann sampler, the rejection probability of z is 1 — p B (z\ {z, z}). It is equal to the posterior likelihood 
p B (z\ {z, z}) of z and to the acceptance probability of z constructed from 5 using the proposal probability q a ({z, z} \z). 
Posterior likelihoods relate to marginal and conditional probabilities via Bayes theorem 

B( r ~t\ |g({M} \zc)P(z c ) . . 

P {Zc\{z,z}) = i; 77 ~ 1X (A2) 



where z c is either z or z. As a result, the Boltzmann sampler leaves invariant the probability distribution $e ({z, z}), 
since it satisfies the detailed balance equation (z c 6 {z, 2} n {z', z'}) 



q({z', z'} \z c )p B (z c \ {z, z})$ e ({z, z}) = q({z, z} \z c )p B (z c \ {z , z'})$ e ({z', z}). 



(A3) 



In Eq. ([A3]), both p B (z c \{z\ z' }) and p B (z c \{z' , z'}) correspond to either an acceptance or a rejection probability. 

Let r Q denote the space of webs sharing a common state at t = ar. The total web space is Tq U Ti and we have 
equipartition of webs in Tq and T\. The free energy differences can be cast in the formii (a 6 {0, 1}) 



3a ({z, z}\z)P e {z)e 0(0 ~ a)w{z) + q a ({z, z}\S)Pg(S)e^ e - a)W{i) VzV5 
p B (z\{z, ~ z})e P(o- a )w iz) + p B ~ z}) e W- a) W { ^ M{z ^ ~ })VzV ~ z 

l 



J(a-8)W(z) , B(a-0)W{Z) 



<& e {{z,z})VzVz 



Since the algorithm sample the distribution $g, the free energy difference gg — g a can be estimated from the logarithm 
of the following estimator 



K n (eW-°» w ) = [e^ Q - 



a /3(a-9)W(2 2fe4 



fe=l 



(A4) 



In Eq. (|A4I) . trial paths z 2 k and z 2 k+i are generated forward and backward, respectively. Hence, we have 
(z2k+a, Z2k+a) € T a . For more details, see Ref. I 111 
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